
tab.path <- "results/table/"
fig.path <- "results/graph/"
fig.name <-  "app_fig_g01_a.pdf"
fig.name2 <-  "app_fig_g01_b.pdf"

#FOR ROBUSTNESS 0.2 to 1 deg cluster x 0.2 to 1 deg cluster
startlat <- min(Re$latgroup)
startlong <- min(Re$longgroup)

NKlus <- c()
Re$y <- 0
Re$x <- 0

for (i in 1:10) {
j <- 0
  while (startlat + i*(j+1) <= max(Re$latgroup)){

   Re$y[Re$latgroup >= startlat + (i)*j 
    & Re$latgroup < startlat  + i*(j+1) ] <- j

  Re$x[Re$longgroup >= startlong + (i)*j 
    & Re$longgroup < startlong  + i*(j+1) ] <- j
  j <- j +1
  }
Re[[paste0("gridklust_",i)]] <- paste0("klus_",Re$y, "_", Re$x)
Nk <- length(unique(Re[[paste0("gridklust_",i)]]))
print(paste("total number of clusters for iteration", i, ":"))
print(Nk)
NKlus <- c(NKlus, Nk)
Re$y <- 0 #reset
Re$x <- 0
}

Kap <- paste(paste("Grid", 1:10/10 , "deg x",  1:10/10 , "deg: "), 
    NKlus)
Kap2 <- paste(Kap, collapse = "\n")


yvars <- c("par.onl.sh")

fcontrols <- paste0("dpost:",controls)

coefs <- 0
se <- 0
klus <- 0
spec <- 0
Nk <- 0

df <- data.frame(klus, coefs, se, spec, Nk)

for ( i in 1:10){
tab.name <- paste0("ols_did_cluster_", i, ".tex")
Klus <- paste0("gridklust_", i)
vars <- c( "factor(dpost)*factor(treat_10)","factor(dpost)" ,"factor(treat_10)")
outvars <- c( "factor.dpost.TRUE.factor.treat_10.TRUE","^factor.dpost.TRUE$" ,"^factor.treat_10.TRUE$")
labvars <- c("DPost X March","DPost", "March")


m1  <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars) ,
      FE = "county" , IV="0", clust = Klus )),
           data= Re)
m2 <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = Klus)),
           data= Re)
m3 <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = Klus)),
           data= subset(Re, year %in% c(1912, 1914)))
m4 <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = Klus)),
           data= subset(Re, year %in% c(1911,1912, 1914)))
m5 <- felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = Klus)),
       data= subset(Re, year %in% c(1911,1912, 1913, 1914) 
        & (pop < 15000 & (dist.roads <= 2 | treat_10 == 1))))

getwd()

regstoadd <- list(m1,m2, m3,m4, m5 )

MSd <- function(variable){
  x1 <- round(mean(variable, na.rm = TRUE), 2)
  x2<-  round(sd(variable, na.rm = TRUE), 2)
  return(c(x1,x2))
}

Stats <- rbind(MSd(Re[,yvars[1]]),
  MSd(Re[,yvars[1]]) ,
  MSd(Re[,yvars[1]]) ,
  MSd(Re[Re$year %in% c(1912, 1914),yvars[1]]),
  MSd(Re[Re$year %in% c(1911, 1912, 1914),yvars[1]]),
  MSd(Re[Re$year %in% c(1911, 1912, 1914) & Re$pop < 10000 ,yvars[1]])
)
####OUT
add.lines <- list(AddLines( 7, "County FE", rep(LMc("Yes"),5)),
                  AddLines( 7, "Incl. 1913", c(rep(LMc("Yes"),2), c(LMc("No"), LMc("No"), LMc("Yes")))),
                  AddLines( 7, "Incl. 1911", c(rep(LMc("Yes"),2), c(LMc("No"), LMc("Yes"), LMc("Yes")))),
                  AddLines( 7, "Controls", c(LMc("No"),  rep(LMc("Yes"),4))),
                  AddLines( 7, "Pop under 15k", c(rep(LMc("No"),4), LMc("Yes"))),
                  AddLines( 7, "Within 2~km of roads ", c(rep(LMc("No"),4), LMc("Yes"))),
                  AddLines( 7, "Mean dep. var." , Stats[,1] ),
                  AddLines( 7, "Sd dep. var." , Stats[,2] )
                 )

fileConn<-file(paste0(tab.path,tab.name))
print(fileConn)
writeLines(stargazer(regstoadd,
                     float = FALSE ,
                     keep= outvars,
                     order = outvars ,
                     covariate.labels= labvars ,
                     star.cutoffs = c(0.1099, 0.0599, 0.0199) ,
                     digits.extra = 0,
                     multicolumn=F,#
                     #column.sep.width="10pt",
                     dep.var.caption = paste("Share of Local Electors" ),
                     dep.var.labels.include = F,
                     font.size = "footnotesize",
                     omit.table.layout = "n",
                     align =T,
                     add.lines = add.lines,
                     digits = 3,
                     intercept.top = T,
                     intercept.bottom = F,
                     keep.stat = c("rsq" , "n") ,
                     omit.stat = c("res.dev","ser") ),
           fileConn)
close(fileConn)


# Make graphs

kall <- "factor(dpost)TRUE:factor(treat_10)TRUE"

coefs <- c(coefficients(m1)[kall],
    coefficients(m2)[kall],
    coefficients(m3)[kall],
    coefficients(m4)[kall],
    coefficients(m5)[kall]
)

se <- c(m1$cse[kall],
    m2$cse[kall],
    m3$cse[kall],
    m4$cse[kall],
    m5$cse[kall]
    )

spec <- 1:5

Nk <- length(unique(Re[[paste0("gridklust_",i)]]))

klus <- rep(i/10,5)

df <- full_join(df, data.frame(klus, coefs, se, spec, Nk))

}
df <- subset(df, klus != 0)

df$fspec <- cut(df$spec, 5, 
    labels = c("Spec (1) - No Controls",
        "Spec (2) - All Sample, all Controls",
        "Spec (3) - Only 1912 and 1914, all Controls", 
        "Spec (4) - Exclude 1913, all Controls", 
        "Spec (5) - Exclude Urban Centers and Parishes far from Roads"
    )
    )

df$up10 <- df$coefs + 1.65*df$se
df$lo10 <- df$coefs - 1.65*df$se
df$up5 <-  df$coefs + 1.96*df$se
df$lo5 <-  df$coefs - 1.96*df$se

kols <- ghibli_palette(name = "MononokeLight")[1:5]
 
p.t <- ggplot(data = df, aes(x = klus)) +
    geom_pointrange(aes( y= coefs, ymin=lo5, ymax = up5, shape = fspec , col = fspec) , 
    size = 1, fatten = 2,  
    position = position_dodge(width = 0.07)) +
    geom_linerange(aes(ymin=lo10, ymax = up10, col = fspec) , 
    size = 3,  alpha = 0.25, 
    position = position_dodge(width = 0.07)) +
    geom_vline(xintercept=0,linetype="dotted", size=0.2) +
    geom_hline(yintercept=0,linetype="dotted", size=0.2) +
    scale_color_manual(values = kols, name = "Specification") +
    scale_shape_manual(values = c(0,1,2,8,16), name = "Specification") +
    scale_x_continuous(breaks=1:10/10) +
    theme_light(base_size = 16) + 
    theme(legend.position = "bottom",
        legend.direction = "vertical") +
    labs(x = "Size of the cluster grid in degrees",
        y = "Effect Post x March")
ggsave(plot = p.t , paste0(fig.path, fig.name), width = 8, height = 6)

p.t2 <- ggplot(data = df, aes(x = klus)) +
    geom_point(aes(y = Nk) , col =  ghibli_palette(name = "MononokeLight")[1]) +
    geom_line(aes(y = Nk) , col =  ghibli_palette(name = "MononokeLight")[1]) +
    scale_x_continuous(breaks=1:10/10) +
    theme_light(base_size = 16) +
   labs(y = "Number of Clusters",
        x = "Size of the cluster grid in degrees")
ggsave(plot = p.t2 , paste0(fig.path, fig.name2), width = 8, height = 2.5)

rm(list=
    setdiff(ls()[sapply(ls(), function(x) any(class(get(x)) == 'data.frame'))],
    c("Re","Ind")))
